Dynamical Selection in Emergent Fermionic Pairing 
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We consider evolution of a Fermi gas in the presence of a time-dependent BCS interaction. The 
pairing amplitude in the emergent BCS state is found to be an oscillatory function of time with 
predictable characteristics. The interplay of linear instability of the unpaired state and nonlinear 
interactions selects periodic soliton trains of a specific form, described by the Jacobi elliptic function 
dn. While the parameters of the soliton train, such as the period, amplitude, and time lag, fluctuate 
among different realizations, the elliptic function form remains robust. The parameter variation is 
accounted for by the fluctuations of particle distribution in the initial unpaired state. 



Nonequilibrium effects in superconducting systems are 
usually described using the notion of local equilibrium 
and quasiparticle distribution [1,2], embodied in the 
time-dependent Ginzburg-Landau equation [3] or the ki- 
netic equation [4]. The recently studied problem of fast 
time dynamics in emergent fermionic pairing [5] describ- 
ing the onset of BCS pairing triggered by an abrupt 
change of interaction lies somewhat outside this classi- 
fication. The new aspect of this problem is the absence 
of a meaningful notion of quasiparticle spectrum. In- 
stead, individual Cooper pair states evolve in a coherent, 
collisionless fashion [6-9] independently of the order pa- 
rameter. The ensuing interesting many-body evolution 
reflects the system inability to transit adiabatically be- 
tween the unpaired and paired states due to their van- 
ishing overlap. The oscillatory mode [5] can be viewed 
as Bloch precession of the pair states in an effective mag- 
netic field defined in terms of the pairing function. 

Tunable time-dependent pairing interaction can be re- 
alized experimentally in ultracold Fermi gases [10-13]. 
Fermionic pairing and superfluidity in these systems have 
been demonstrated recently using magnetically tunable 
Feshbach resonances [14-16]. In these systems, by vary- 
ing the magnetic field detuning from resonance, both the 
strength and the sign of inter-particle interaction can be 
changed on a time scale shorter than the intrinsic times 
of fermions, set by Ep or the elastic collision rate. In 
Refs. [18,20,21] the picture of time-dependent pairing [5] 
was extended to the problem of resonant atom-molecule 
coupling near a Feshbach resonance. 

Another interesting class of tunable systems are the hy- 
brid semiconductor-superconductor (S/Sm/S) and SNS 
structures [24-28] with proximity-induced Josephson ef- 
fect. In the S/Sm/S systems [24,25], the superfluid den- 
sity and Josephson critical current can be tuned by elec- 
tric field applied on the gates to the semiconductor. Sim- 
ilarly, the SNS structures [26-28] are tunable by current 
applied to the normal region. 

Superconductor dynamics is controlled by several ki- 
netic time scales [2]. The most important for us is the 
collisionless regime [6-9] characterized by the time 

r A = ft/ A (1) 

with A the equilibrium BCS gap. The time ta manifests 



itself in the BCS instability growth rate [7-9] at T away 
from T c , as well as in the frequency of pairing amplitude 
oscillations [6]. Another time scale is set by the Lan- 
dau Fermi liquid elastic collision time, r o T 1 oc max[e 2 , T 2 ] 
with e the energy relative to Ep. For typical energy 
e ~ A, in the weak coupling regime A -C Ep, we es- 
timate r e i 3> ta. For slowly varying interaction, the 
kinetic equation [4], or the time-dependent Ginzburg- 
Landau equation [3] can be employed. In contrast, here 
we are interested in the situation when the interaction is 
turned on abruptly on the scale ta , so that the hierarchy 
of times is 



T < T A < Tel, 



(2) 



where To describes the interaction time dependence. In 
this case, the time interval < t < t c \ is controlled by 
coherent pair dynamics, while energy relaxation can be 
ignored. This regime can be described by the truncated 
BCS Hamiltonian, which accounts for coherent Cooper 
pair transitions (p, — p) — * (p', — p'), ignoring other, 
more slow, processes of quasiparticle scattering. In this 
article we study the resulting collisionless evolution, us- 
ing a combination of analytical and numerical methods. 

The time domain picture of BCS pairing buildup, pro- 
posed in Ref. [5], involves dynamics of individual pair 
states, selfconsistently coupled to the pairing amplitude. 
The latter exhibits unharmonic oscillations, undamped 
at t < t c i. While several recent publications agree in 
general with this scenario [19-22], some of the issues re- 
main unsettled. Ref. [19] describes a simulation of the 
collisionless BCS dynamics exhibiting a damped oscilla- 
tory time dependence, with the damping of unspecified 
origin. These results, as well as those of Ref. [23], are 
clearly different from the behavior found in Ref. [5] . The 
most probable origin of the discrepancy, in our view, is 
due to the numerical method [19]. 

Below we provide extensive numerical evidence which 
confirms the conclusions of Ref. [5]. We show that the 
mean field solutions, the soliton trains described by el- 
liptic dn function, are singled out by the BCS evolution, 
and appear for generic initial conditions in a stable and 
robust way. Estimates of noise made in Sec. IV show that 
the mean field theory approach holds when the level spac- 
ing in a relevant volume is small enough. Our analysis 
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indicates that the moan field results describe the dynam- 
ics of |A| in a wide temperature range not only in a finite, 
but also, locally, in an infinite system. 

The physical reason for a specific solution to be se- 
lected by pairing dynamics is due to the properties of the 
BCS instability. Linearization over the unpaired state of 
a Fermi gas in a finite system, discussed in Sec. I, shows 
that, modulo phase degeneracy, there is just one unsta- 
ble normal mode corresponding to perturbation growth, 
while other modes do not grow. As a result, although 
the initial state is perturbed by thermal fluctuations in a 
completely random fashion, the BCS dynamics amplifies 
only one specific perturbation leading to an oscillatory 
time dependence with predictable characteristics. 



BCS INSTABILITY OF THE UNPAIRED 
FERMI GAS 



A(t) = A^> p (W)- 



(6) 



The amplitudes u p (t), v p (t) can be obtained from the 
Bogoliubov-deGcnnes equation 



A 



A* 



(7) 



solved together with the selfconsistency condition (6). 

We recall that the unpaired state is a selfconsistent, al- 
beit unstable, solution of Eqs. (7), (6) with A = 0, T = 0: 



with 



uZ>(t) = e-^(e p ), v p °\t) = e ^p e -p*0(-e p ), (8) 



mdom phase. The stability analysis [8] 
shows that the deviation from the unpaired state grows 
as A(t) oc e 7 *e~^*, with linearized amplitudes 



The evolution of the Fermi gas with time- varying pair- 
ing coupling can be described by the BCS Hamiltonian 



a p,<j a P,<? 



2 



2^ a PT a - 



-ql a qT> ( 3 ) 



where a 



are the canonical fermion operators, and 



a =t,4 is generalized spin. The time dependence of A, 
as well as the resulting time dependence of the system 
state, is assumed to be fast on the scale of quasiparticle 
elastic collisions and energy relaxation, r e i, allowing us 
to ignore the latter and consider the coherent dynamics 
defined by (3). The simplest time dependence which we 
shall be most interested in, is described by the coupling 
turned on abruptly, from X(t < 0) = to X(t > 0) = A. 

The interaction switching, while abrupt and nonadi- 
abatic, must also be gentle enough not to overheat the 
Fermi system. The analysis of energy production due to 
two-particle scattering in the presence of time-dependent 
coupling, described in Appendix, obtains an estimate for 
the effective temperature 



j.2 



off 



(4) 



The 'no overheat' condition T e g <C A can thus be stated 
as EpTo 3> (An/Ao) 2 / 3 , which is compatible with the 
nonadiabaticity requirement to < ta. 

Our treatment of the problem (3) will focus on the 
time-dependent generalization of the BCS state [5] 



*(*)) =U(u P (t)+v p (t)a+^a: 



p. I 



|0>. 



(5) 



The Bogoliubov mean field approach, which gives a state 
of the form (5), relies on the 'infinite range' form of 
the pairing interaction in (3) owing to equal coupling 
strength of all (p, — p), (q, — q). Since the latter docs 
not depend on the system being in equilibrium, one can 
introduce a time-dependent mean field pairing function 



5u p (t) = - 



ry - 2e p + u ' 



j. j(t) - A-(*)4 0) (*) (9) 
pl ) " i 1 + 2e p -u W 



The growth exponent 7 and the constant u, combined 
into a complex number ( = to + vy, are defined by the 
selfconsistency condition of the linearized problem: 



sgne p 
2e P -C' 



(10) 



This equation has a pair of complex conjugate solutions 
£, (*■ From the similarity between Eq. (10) and the BCS 
gap equation one expects the exponent 7 to be close to 
the BCS gap value A in equilibrium, in which case the 
time constant ta — 7 _1 for initial pairing buildup is of 
the order of Ag" 1 . 

It will be convenient for us to introduce here the con- 
stant density of states approximation, 



f(Cp) 



0, 



else 



<\W 



(11) 



used throughout this article. In our numerical study we 
use a total N » 1 equally spaced discrete states dis- 
tributed evenly in a finite size band, 



-W/2 < e p < W/2, 



(12) 



with the level spacing Se = W/N. Although somewhat 
artificial, the model (12) introduces a convenient simpli- 
fication due to the particle-hole symmetry. In this case, 
the chemical potential is locked at e = independent of 
the interaction strength. As a result, the instability prob- 
lem (18) possesses an eigenvalue with a purely imaginary 
( = fy, satisfying 



2|e r 



4e p + 7 2 ' 



(13) 
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and u> = 0. Interestingly, if the bandwidth W is much 
larger than the BCS gap Ao at T = 0, the value of the ex- 
ponent 7 obtained from (10) coincides with Ao. Indeed, 
Eq.(13) in this case gives 



1 = 2Ai/ 



W ' 2 2ede 1 W 2 + 7 2 M A , 



This equation is reminiscent of the BCS gap equation 



1 



1 



W/2 



de 



W/2 y/e 2 + A 2 



= Xuq sinh 



W 

2A7 



becoming identical to it in the weak coupling limit of 
fermion energy band wide compared to the gap, W ^> 
7, Ao [cf. Eq.(41) below]. The condition u> = is vi- 
olated in the absence of particle-hole symmetry. The 
equality 7 = Ao holds only at T = 0. 

Let us briefly discuss how the analysis is modified in 
the case of a discrete spectrum (12). With the sum in 
Eq.(10) running over N levels, we obtain a polynomial 
of order N, having total N roots. Simple analysis shows 
that N — 2 roots are real and the remaining two are a 
complex conjugate pair £, with values close to that 
obtained for continuous spectrum. Accordingly, only the 
normal modes corresponding to £, (* are relevant for the 
instability while the other N — 2 modes correspond to 
the perturbations which do not grow and remain small. 
This conclusion is valid only at the times described by 
the linearized BCS dynamics. The situation at longer 
times, which is more delicate, will be discussed below. 



where w p , c 
occupancy: 



"p«p,T + c P a ^ P ,l + v P a tl a -pA ) 1°)' ( 15 ) 



equal zero or one depending on the 
,c p ,c' p ,v p ) = (1,0,0,0), ...,(0,0,0,1). 



The average values are given by occupation probabilities: 



''Pi 



f =(l-«p) 2 



"Pi 



12 -n 2 p , 



J2 = | fj |2 = 



; -P 



= n p (l - n F 



where n p — l/(c p + 1) is the Fermi function, and the 
quantities n p , n p (l — n p ) and (1 — n p ) 2 describe double, 
single and zero occupancy probability of the two-fermion 
states (p, — p) in the unpaired system. As we argue be- 
low, while the product state (15), written as a finite tem- 
perature generalization of (5), is not the most general 
fermion state, it is adequate for our problem. 

The product state (15) is suitable for simulation, since 
Bogoliubov-deGennes dynamics (7) couples only u p and 
v p independently for each p, preserving the product 
form. At the same time, the parts of (15) with single 
occupancy are decoupled from the collisionlcss pair dy- 
namics (3). Indeed, the Hamiltonian (3) gives zero when 
applied to singly occupied pair states aj||0), a+ k JO), 
irrespective of the occupancy of other pair states. One 
can thus identify a subspace in the full Hilbert space, 
which is spanned by all combinations of pair states of 
occupancy zero and two. The latter have the form 



)gcncral — a k,| fl -k,l I ^) ' 



(16) 




FIG. 1. Temperature dependence of the BCS instability 
growth rate 7 as obtained from Eq.(18) for the constant den- 
sity of states model (11), (12) with the coupling A such that 
Ao/W = 1/5. Note that 7 coincides with the BCS gap Ao at 
T = 0, up to a correction small as (j/W) 2 at large bandwidth 
W. Near T c — ire~ c Ao — 1.764Ao, the exponent 7 vanishes 
linearly in T c - T. 

Let us now consider the instability at finite tempera- 
ture. The initial state describing Fermi gas at T > can 
be tentatively chosen as 



with the product taken over the states (k, — k) of occu- 
pancy two whereby the states of occupancy one are ex- 
cluded from the vacuum |0). The states (16) are mapped 
by (3) onto the states of a similar form, thereby defining 
a full representation of the Hamiltonian. 

Fortunately, one can bypass the combinatorics of (16) 
and simplify the state by employing Bogoliubov mean 
field approach which is exact for the Hamiltonian (3). In 
the mean field framework, the general state is replaced 
by a product state. Indeed, since the pairing amplitude 
A(t), describing cumulative effect of all pairs, is a c- 
number, the dynamics (7) does not generate correlations 
between different pair states. At the same time, any cor- 
relations present at t = are dephased by the dynamics 
itself, described by time-dependent 2x2 evolution ma- 
trices, different for each pair (p, — p). This argument, 
which will be refined in Sec.V, allows to replace the gen- 
eral states (16) by a simpler state of the product form 
(15) with c p = c' p = 0. On a mean field level, the corre- 
lations between different (p, — p), (k, — k) do not matter. 

The effect of Pauli blocking which eliminates the states 
of occupancy one can be taken into account by going back 
to the T = state (5) with u p and v p chosen as 



(*) 



ept (l-« P ), 



'P' 



(17) 
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with random phase <f> p . The reduced norm |u p | 2 + |f p | = 
n p + (l— rip) 2 < 1 reflects that at finite temperature some 
pair states, populated by just one particle, are decoupled 
from the dynamics. Near the Fermi level, at |e p | <C T, the 
two-particle states with double or zero occupancy have 
the probability 1/4 each, so that |w p | 2 + \v p \ 2 = 1/2. 
Outside this interval, |e p | > T, the blocking is practi- 
cally absent and the norm approaches one. Below we 
shall use the state (5) with u p , v p given by (17) as initial 
condition for the simulation. 

The choice (17), while somewhat ad hoc, has an addi- 
tional advantage over choosing u p = 1, 0, v p = 0, 1. With 
both u p and v p nonzero with random relative phase, the 
state (5) exhibits pairing fluctuations, providing a "seed" 
for the BCS instability in the simulation. The results of 
the latter indicate that this choice of the initial state is 
general enough for the instability to fully play out. As we 
shall see in Sec.V, where a spin 1/2 formulation of BCS 
dynamics is discussed, u p and v p can be understood as a 
two-component spinor. This means that the state of the 
entire system can indeed be chosen in the product form, 
with the initial values u p , v p having random modulus, 
not just phase. This difference, however, is inessential, 
since the modulus of u pi v p is quickly randomized by the 
dynamics itself. 

Returning to the analysis of the instability, lineariza- 
tion of (7) over the finite temperature state (17) obtains 
a time-dependent perturbation of the form Eq.(9), and a 
generalization of Eq.(10), with sgne p replaced by 

|u p 0) | 2 - |4 0) | 2 = 1 - 2n p = tanhi/3e p . 
The resulting equation, 

1 = A E^T' A(t)«e-«*, (18) 

has nonzero solutions £, £* below the BCS transition, 
T < T c , which are purely imaginary, £, £* = in the 
case of particle-hole symmetry (11), (12). The growth 
exponent 7 vanishes as T — > T c , as shown in Fig.l. This 
can be verified by noting that Eq.(18) yields infinitcs- 
imally small 7 at T = T c , since at £ = it becomes 
identical to the BCS equation for T c . 

Having established the form of the unstable mode (9) 
obtained from linearization, let us now estimate the time 
range for which this analysis is accurate. The denomi- 
nator in Eq.(9) is larger than the numerator as long as 
A(t) < 7, with 7 w A(T), the BCS gap. The initial 
value t] = A(t = 0), nonzero due to fluctuations in the 
unpaired Fermi system, is small in 1/N: 

V = \J2u P v* p = \Y / n P (l-n p )e i *». (19) 
p p 

The sum (19) is controlled by about T/Se terms with 
l e pl ^ T, uncorrelated in phase. From the central limit 
theorem argument, we estimate 77 ~ \v\fTbk by order 



of magnitude. The condition A(t) < 7, with exponen- 
tially growing A(t) = 7/e 7 *, defines the time interval 
< t < 7 _1 0(7/77) in which the evolution is described 
by the linearized problem. 

Our aim will be to gain insight in the behavior at 
later times. We rely on the nondissipative character 
of the Bogoliubov-deGennes dynamics (7), (6), manifest, 
for instance, in the energy E = {t)\TL\^ {t)) conserva- 
tion throughout the evolution. Since for the unpaired 
state the energy exceeds its value in the ground state 
by the BCS condensation energy, the equilibrium cannot 
be reached without collisions. Eqs. (7), (6) hold at times 
shorter than quasiparticle thcrmalization time r e i. For 
temperatures away from T c , the time t c \, evaluated at 
e ~ A , is long compared to the instability growth time, 

r e i > 7 _1 - (20) 

This means that the linear instability phase is followed 
by a long collisionless nonlinear phase. Below we show 
that the evolution governed by Eqs. (7), (6) is described 
by soliton-likc pulses in A(t). We obtain a family of exact 
solutions of the form of single solitons and soliton trains, 
and compare it with simulations. 

We briefly note that the importance of coherent dy- 
namics of individual Cooper pairs (5) in the time evolu- 
tion of a paired state has been understood a long time ago 
in the context of the discussion of the validity of the time- 
dependent Ginzburg-Landau equation approach [7,8,3]. 
It has been pointed out by Gorkov and Eliashbcrg [3] that 
the time-dependent pairing function A is generally insuf- 
ficient to describe the evolution. For such a description 
to be consistent, the pair breaking and energy relaxation 
must be fast compared to the time scale ta of change of 
A. This can be realized only close enough to the transi- 
tion, or in superconductors with magnetic impurities [3] , 
where the inequality (20) is violated. Except these spe- 
cial situations, however, the Cooper pairs are not slaved 
to the time-dependent A(i), and the dynamics of each 
pair has to be treated individually, via Eq. (7). 

II. BOGOLIUBOV-DEGENNES EQUATION AS A 
BLOCH EQUATION 

To analyze the nonlinear BCS dynamics we reformu- 
late the Bogoliubov approach, bringing it to a form 
amenable to analytic and numerical treatment. We show 
that the evolution of time-dependent amplitudes u p (t), 
v p (t) , governed by the Bogoliubov-deGennes equation (7) 
with the selfconsistency condition (6), can be cast in the 
form of a Bloch equation for auxiliary variables. This is 
achieved by by introducing a new set of variables, 

fp = 2u p v;, g P = \u p \ 2 -\v;\ 2 . (21) 

Applied to the quantities (21), Eq.(7) gives a system of 
coupled equations: 
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^ = -2ie P / P + 2tAg p , -|E = iA*f p - iA%. (22) 

These are nothing but the Gorkov equations [29] for the 
Green's functions G and F. In the form (22), using 
momentum-dependent quantities f p , g pi these equations 
were first written by Volkov and Kogan [6]. 

The dynamics (22) is to be supplemented by the 
Gorkov sclfconsistency relation [29], 

A=^£ip> (23) 
p 

which defines A(i) through the values of evolving f p (t). 
The initial conditions 

4 W*p(l-np)np, 5 W = l-2n p , (24) 

correspond to the unpaired Fermi gas state (17). 

To gain insight in the behavior of Eqs.(22), it will be 
convenient to introduce the Bloch representation 

r p = (ri,r 2 ,r 3 )p, r 1 +ir 2 =f p , r 3 = g p . (25) 

The norm of r p is given by 

|r P | = y/\f P \ 2 +9 v = 4 + 0- "p) 2 - ( 26 ) 

Remarkably, after rewriting Eq.(22) in terms of r p , it 
assumes the form of a Bloch equation: 

dr 

= 2b P x r p , (27) 

where the "magnetic field" b p = —(A', A", e p ) has time- 
dependent x and y component satisfying the selfconsis- 
tency condition (23). The norm of the Bloch vectors r p , 
given by Eq.(26), is less than one, which describes the 
effect of Pauli blocking, i.e. decoupling of the states with 
single occupancy from the BCS dynamics (see Sec. I). 




0.2 0.4 0.6 0.8 1 



Temperature T/T 

FIG. 2. Temperature dependence of the pairing amplitude 
A for the stationary state (28), (29) obtained from the un- 
paired state by adiabatic increase of coupling. The equilib- 
rium BCS gap is shown for comparison. 



Before exploring the dynamics, let us inspect the sta- 
tionary states of Eq.(27). The latter are described by 
the Bloch spins aligned with the magnetic field axis, 
l p = — bp/|b p |. In this case we have 

(h+ih) P = -F A Z 3 , p = r-^-^ , (28) 

^ p +i A i 2 ^Ap+i A i 2 

and the selfconsistency condition (23), which determines 
the stationary value of A, takes the form 

i = A tanh^/3|ep| ^ 

2 V^| + |aF 

The numerator, tanh(|/3|e p |) = |1 — 2n p |, is the length 
of the Bloch vector f p averaged over a group of levels 
with nearly equal energies, |e Q — ep\ <C A. The averag- 
ing, applied to the initial r p values (24), eliminates the 
transverse part of (24), containing random phases 9 P , 
while leaving the longitudinal part intact. The Bloch dy- 
namics is unitary with respect to each Bloch vector r p 
and, in particular, is linear and preserves the norm. As 
a result, the averaged vectors r p will evolve according to 
the same Bloch equations, albeit having a smaller norm 
|f p | = |l-2np|<|r p |. 

The equation (29) is different from the BCS gap equa- 
tion which contains tanh(i/3(e 2 + A 2 ) 1 / 2 ) instead of 
tanh(i/?|e p |), except T = 0, when these equations co- 
incide since tanh = ±1 in both cases. Thus at temper- 
atures < T < T c , as Fig. 2 illustrates, Eq.(29) predicts 
the stationary value of A below the BCS gap scale. The 
temperature at which A vanishes coincides with the BCS 
critical temperature, since the condition (29) at A = is 
identical to the BCS equation for T c . 

To clarify the character of the states (28), (29), one can 
make the following observations. The only difference here 
from the BCS theory is due to incomplete equilibrium, 
owing to the singly occupied states being Pauli-blocked 
from the dynamics controlled by (3). Indeed, the trun- 
cated BCS Hamiltonian (3) accounts only for the colli- 
sionless pair dynamics, but not for single particle scatter- 
ing and relaxation. The latter processes have character- 
istic rates set by the two-particle collisions, l/r e \. Since 
r e i is larger than ta = 7 , the approach accounting only 
for the coherent dynamics, but not for relaxation, is valid 
in a relatively large time interval < t < T e \. 

The stationary nonequilibrium states (28), (29) can be 
realized when the coupling constant A increases as a func- 
tion of time slowly on the scale of ta, i.e. the condition 
(2) is replaced by 

TA < T < T C 1. 

In this case, each Bloch vector r p follows the direction of 
the field b p (t), maintaining constant projection on the 
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b p (t) axis equal to 1 — 2n p on average. During the evo- 
lution, the value A(t) is determined by the selfconsis- 
tency condition (23). The amplitudes of the pair states 
with occupancy zero and two thereby become slaved to 
the adiabatic dynamics A(£), evolving according to (28), 
(29). At the same time, the amplitudes with occupancy 
one remain decoupled, and do not evolve at times t < t c \. 
Such a behavior can be seen as a result of the evolution 
which is simultaneously adiabatic in the pair sector, and 
totally nonadiabatic in the single particle sector. 



u> random in sign, constant for each realization, equal 
twice the chemical potential. In addition, we observe 
noise superimposed on the linear time dependence of the 
phase, which we shall discuss below at the end of this 
section. As noted in Ref. [5], finite w can be eliminated 
by a gauge transformation which shifts single particle en- 
ergies by uj/2, e p = e p — lo/2, thereby making A real. In 
the Bloch equation language, this is equivalent to consid- 
ering the problem in a Larmor frame rotating about the 
z axis with frequency u>. Having this in mind, below we 
shall focus on the behavior of the modulus |A|. 



III. OSCILLATORY SOLUTIONS, ANALYTICAL 
AND NUMERICAL 



Here we consider the nonlinear BCS dynamics de- 
scribed by the Bloch equation (27) and the selfconsis- 
tency relation (23) at the times after the instability sets 
in. Eq. (27) is quite easy to simulate, since it is linear in 
r p and is written for classical, rather than quantum vari- 
ables. The initial state, Eq.(24), describing free fermions 
at a finite temperature, is 



i<f>p 



n, p +tr2,p 



cosh(i/3e p ) 



r 3 , p = tanh(~/3e p ), (30) 



with uncorrelated phases, uniformly distributed in the 



interval < 



< 2-7T. This form of initial conditions 



corresponds to the amplitudes u p , v p of the form (17). 
To avoid confusion with temporal characteristics, such as 
the oscillation period T, hereafter we shall use the inverse 
temperature j3, unless explicitly stated otherwise. 

The dynamics (27), (23) was obtained from 3N cou- 
pled differential equations with the initial conditions (30) , 
with N large enough to ensure proximity to the continual 
limit. The numerics was executed using the Runge-Kutta 
method with precision 0(dt 5 ). The time step dt was var- 
ied over a range of values to test numerical accuracy. We 
found that the step dt = 0.01/Ao typically provides suf- 
ficient precision over the time interval of interest. 

As Fig. 3 illustrates, a straightforward simulation gen- 
erates a surprisingly regular, oscillatory time dependence 
A(t) which appears to be a periodic function of time. Af- 
ter initial exponential growth, controlled by the instabil- 
ity discussed in Sec. I, we observe an essentially periodic 
time dependence, characterized by equally spaced peaks 
of identical shape. We can thus define the temporal pe- 
riod T, the time lag r, and the amplitude which we de- 
note by A + for the reasons to become clear below. These 
notations are marked in Fig. 3. 

It might seem that the particle-hole symmetry of the 
density of states (11), (12) would enforce zero chemical 
potential, regardless of BCS interaction strength. In the 
simulation, however, we observe nonzero values of chem- 
ical potential due to particle-hole imbalance caused by 
thermal fluctuations in the initial distribution n p . The 
complex pairing amplitude exhibits a phase growing lin- 
early with time, A(t) = e - ^W|A(t)|, <j)(t) oc ut, with 
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FIG. 3. Time dependence of the pairing amplitude A 
recorded from simulation with N = 10 5 states (11), (12) at 
temperature T = 0.7T C (/3 = 2.5Ao) with the initial condi- 
tions (30). The coupling constant A was chosen to have the 
BCS gap A = W/5. 

Interestingly, the obtained time dependence A(t) can 
be fitted extremely accurately to the analytic solution 
found in Ref. [5] for noiseless initial conditions. The lat- 
ter is given by a Jacobi elliptic function, periodic in time, 

A(t) = A+dn(A + (i-r),fc), k 2 = 1 - A 2 _/A 2 + , (31) 

with A + the amplitude, t the time lag, and A_ the min- 
imal value. We recall that the function u = dn (x, k) is 
obtained by inversion of an elliptic integral: 



dv! 



(32) 



/„ y/{l-u' 2 )(k 2 -l + u' 2 ) 

The function (31) satisfies the differential equation 

(dA/dt) 2 = (A 2 + -A 2 )(A 2 -A 2 _) (33) 

with A± being the extremal values: A_ < A(t) < A + . 
The period of the function (31) is given by the complete 
elliptic integral of the 1st kind: 



2 r /2 d<f> 
A+ Jo y/l - k 2 sin 2 , 



(34) 



G 



For sparse soliton trains, TA + 3> 1, this expression sim- 
plifies toT = ln(4A+/A_). 

The time dependence of the Bloch vectors r p (t) can be 
obtained from the ansatz 



z p = ApA(t) + iBpA(t), r 3 , p = C P A" - D p . (35) 



The role of this equation is similar to the BCS gap equa- 
tion. The only difference is that it fixes one of the two 
constants A±, leaving the other one free. 

The motivation to consider this particular solution can 
be seen from the behavior of the elliptic integral (32) at 
A_ = 0. In this case, we obtain a single soliton 



Eqs.(22) are satisfied by (35) provided A p — 2e p B p and 



Bp = — C p . The normalization condition r\ - 



-4 = i 



thereby turns into Eq.(33), the same for all p, yielding 
the relation between C p , D p and A±: 



D 2 - 1 D 
p - a2 a 2 2— £ 
Cl - A - A +' l n 



4e p + Ai 



(36) 



Here C p and D p , and likewise A p and B p , depend on e p , 
while the quantities A± are the same for all p. 
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FIG. 4. Temperature dependence of the soliton train am- 
plitude A+, obtained from the selfconsistency condition (37) 
at different ratios A_/A+. Note that at A_ = the ampli- 
tude A+ equals the BCS instability growth increment 7 (see 
Fig.l), while at A_ = A+ the result for the stationary state 
is reproduced (see Fig. 2). 

A special property of the ansatz (35) which makes 
it compatible with the selfconsistency condition (23), is 
that / p = A p A(t) has the same time dependence as the 
left hand side of Eq.(23). Therefore, the selfconsistency 
will hold at all times provided that the quantities A± are 
chosen to satisfy 1 = |^ |r p |A p . Here the averaging 
f p over a group of levels with close energies is performed 
in the same way as in the derivation of Eq.(29). (As 
above, the averaging is compatible with unitary evolu- 
tion due the linear character of Bloch dynamics.) After 
substituting the expressions for A p in terms of A±, and 
r p = 1 — 2n p , we obtain 



2e p tanh(i/3e p ) 



A%y 



4A 2 _A 2 h 



1/2- 



(37) 



A(t) 



cosh A+(t — to) 



(38) 



At large negative time, Eq.(38) describes exponential 
growth of A. Furthermore, Eq.(37) at A_ = is identi- 
cal to the condition (18) for the instability growth rate, 
so that A + = 7. Thus the single soliton solution (38) 
describes the nonlinear evolution following the linear in- 
stability regime. The nonmonotonic behavior of A(i), 
first growing and then decreasing to zero, can be under- 
stood as a result of energy mismatch of the BCS ground 
state and the unpaired state: Energy conservation in the 
collisionless dynamics prevents system to evolve to the 
ground state with lower energy. 

Remarkably, while these solutions appear to be very 
special, they are robust in the presence of noise. Below 
we study the instability of Fermi gas at finite tempera- 
ture, and find that the time dependence survives ther- 
mal fluctuations in the initial state. The reason for such 
a behavior is owing to the property of BCS instability, 
discussed in Sec. I, to develop through a single unstable 
mode. As a result, only the fluctuation in the initial state 
along the unstable direction is amplified by BCS dynam- 
ics, while other fluctuations remain small, providing a 
selection mechanism for the solutions (37). 

Returning to the analysis of soliton trains (31), we 
note that the ratio r — A_/A + controls the inter- 
soliton time separation. Different regimes can be qual- 
itatively understood by noting that A varies in the in- 
terval A_ < A(t) < A + . For r increasing from to 
1, the soliton train period T decreases, making the soli- 
tons overlap stronger and gradually merge, turning into 
weak harmonic oscillations with frequency 2A + as A_ 
approaches A + (see Fig. 1 of Ref. [5]). 

As a function of temperature, the quantity A + varies 
from the value close to the BCS gap Ao at T = 0, to zero 
at T = T c (see Fig.4). At A_ < A+, Eq. (37) turns into 
Eq. (10) which, as we found above, defines the amplitude 
of a single soliton (38). In the opposite limit, A_ — > A + , 
Eq. (37) coincides with Eq.(29) for the stationary state. 

To understand the behavior of A + in more detail, let 
us analyze the selfconsistency condition (37) for the sym- 
metric band of states (12). We first consider zero tem- 
perature, when tanh i/3e = sgne. The integral (37), eval- 
uated using variable substitution x — 4e 2 , gives 



dx 



(x + A\ + A 2 _) 2 -4A2_A 



cosh 



_ x w 2 + (r 2 + 1)A 



2rA^ 



mi. 

r 



(39) 
(40) 
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Substituting this in Eq. (37) and solving it, we obtain 



+ 2(e 1 /A_ r 2 e -i/A) sinh i ' v ; 

where here and below the density of states is absorbed 
in the coupling, \vq — > A. At r = 1 we recover the BCS 
gap for the symmetric band (12), Ao = M / /2sinh j-. At 
r < 1 we obtain a value somewhat below A , the differ- 
ence being small as (1 — r 2 )e~ 2 / A Ao. This explains small 



departure from Ao seen in Figs. 1,4 at T — 0, as well as 
its absence in Fig. 2. 

The simulated dynamics A(i) appears to be periodic 
(or very close to it) over the entire time interval of the 
simulation. A nearly perfect fit to A(i), as illustrated 
in Figs. 5, 6, is provided by the elliptic function (31) with 
the same period and amplitude. The numerical and ana- 
lytic functions are found to agree to accuracy better than 
1(T 4 for (3 = 100/Ao (Fig. 5 inset, top panel) and 1CT 3 
for j3 = 10/Aq (Fig. 6 inset, top panel). 
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FIG. 5. Top panel: Comparison of the time dependence A(t) obtained from BCS/Bloch dynamics (27), (23) for N = 10 3 
spins at temperature T = 10 _2 Ao (blue curve) to the analytic soliton train solution (31) of the same amplitude and period 
(green curve). The difference of the simulated and analytic A(t) is shown in red. (The initial conditions (30) and parameters W, 
Ao are the same as in Fig. 3.) Lower panels: The pair distributions of the soliton train parameters for 500 different realizations: 
the time lag and period (left); the period and amplitude (right). 



While each realization A(i) is essentially a perfectly 
periodic function of time of the form (31), the parameters 
such as the period T, the time lag r, and the amplitude 
A_|_ exhibit significant variations from one realization to 
another. To explore this phenomenon, we generated a 
large number (500) of different realizations, and for each 
of them determined the values T, r and A + from fitting 
to the elliptic function (31). Figs. 5, 6 display the result- 



ing pair distributions as a set of points in the (T, r) and 
(T, A_|_) planes, one point per realization. 

These results lead to a number of interesting obser- 
vations. First, as one would expect, the distribution is 
less noisy at lower temperature (Fig. 5). Second, while at 
j3 = 10/ Ao the points are scattered over a 2d region, at 
(3 = 100/Ao each distribution collapses on a Id curve, 
indicating a specific relation between T, r and A + at the 
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lower temperature. 

The pair distribution of the period and the time lag 
tends to cluster around the straight line 

t = T/2. (42) 

This can be understood from the linear stability analysis 
of Sec. I. Indeed, there we found that the linearized BCS 
problem has only two eigenvalues outside the unit circle, 
£ = lo + 17 and = lu — 17. The projections of the initial 
unpaired state on these two vectors are close by order of 
magnitude. Now, let us take into account that the ellip- 
tic function (31) at large period TA + ^> 1 represents a 
train of well-separated 1/cosh solitons (38): 

A(^E coshA A ;_ U +^ T )^ A +=7 , (43) 



t n = nT+T. In between the solitons this function is given 
by a sum of exponential tails of the nearest solitons: 

A(i n < t < t n+1 ) = 7( e -^-*") + e 7(*-*« +1 )). 

The amplitudes of the two terms, taken at some t in the 
interval (i„,i„ + i), should match the C,* projections of 
the initial state for the numerical and analytical solu- 
tion to coincide. Since the terms e -7 '* - *"-', e 7 ^~ t "+ 1 ^ ) 
are equal at the midpoint to = ^(t n + t n+ i), the time t 
defined by the matching condition must be close to to- 
This suggests that the time lag r should indeed be close 
to half a period, with the product (2r — T)*y of order one 
and random in sign. This conclusion is consistent with 
our observations at different temperatures (Figs. 5, 6). 
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FIG. 6. Same as in Fig.5 for higher temperature T = 10 Ao. The simulated time dependence A(t) can be accurately fitted 
to the analytic solution (31), with the distribution of the period, amplitude and time lag somewhat broader than in Fig.5. 



To complete the description of the behavior of A, ob- 
tained in simulation, here we analyze the phase dynam- 
ics. The phase (f> = argA(i) recorded in the simulation 
exhibits an approximately linear time dependence, as il- 
lustrated in Fig. 7. As discussed above, in a system with 



perfect particle-hole symmetry one expects the chemical 
potential to be pinned to the band center, in which case 
the condition d<j)/dt = 2/i = would make the phase 
time-independent. The observed linear behavior can be 
explained by particle-hole imbalance due to fluctuations 
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of particle distribution n p in the initial state. These fluc- 
tuations result in nonzero chemical potential of random 
sign. The fluctuations in /i caused by random occupancy 
will be estimated in Sec. IV, Eq.(50.ii). The magnitude 
and temperature dependence of the fluctuations is found 
to be consistent with observations. 



linear does not lead to noticeable deviation in |A| from 
the elliptic function time dependence. 



IV. NOISE DUE TO OCCUPANCY 
FLUCTUATIONS 
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FIG. 7. Top: The phase of the pairing amplitude versus 
time for the soliton train in Fig. 6. Bottom: Pairing ampli- 
tude A(t) — A x (t) + iA y (t) trajectory in the complex plane. 
The phase arg A(t) is a linear function of time superimposed 
with noise. Each radial line in the left panel corresponds 
to a soliton, marked according to their order in the time se- 
quence. Phase shift between solitons translates into rotation 
by a constant angle. The right panel shows the behavior near 
the origin, allowing to trace the order of different solitons. 



The noise superimposed on the linear time dependence 
4>{t) has several interesting features. First, the fluctua- 
tions about the linear dependence show no sign of phase 
diffusion, since 5(f> does not grow. Instead, they can be 
described as a periodic, or quasiperiodic, arrangement of 
steps connected by kinks. By comparing to Fig. 6 which 
depicts I A I for this simulation, we see that each step is as- 
sociated with a soliton, while the kinks occur in between 
the solitons. This behavior can be partially understood 
by analyzing the trajectory A(t) in a complex plane, dis- 
played in Fig. 7. Each radial straight line in this plot cor- 
responds to a soliton, making it apparent that the phase 
variation occurs mainly in between the solitons. Interest- 
ingly, the deviation of the phase time dependence from 



The robustness of the elliptic function (31) accompa- 
nied by sometimes significant variations of the parame- 
ters T, r and A + among different realizations may seem 
surprising. To get insight into the origin of this behavior 
we consider variations in the initial conditions, which can 
be attributed to fluctuations of the pair states occupancy 
(15) at t = 0. The latter is due to thermodynamic fluctu- 
ations of fermion occupancy n p , and can be interpreted 
more intuitively as temperature fluctuations. 

The first thing we note is that the existence of the 
soliton solutions, as well as their analytic form, is not de- 
pendent upon the details of the energy distribution n p , 
provided the state is unpaired, i.e. there is no coherence 
in the amplitudes u p , v p at different p. The main differ- 
ence arising for the more general distribution is possible 
lack of particle-hole symmetry relative to Ep. In this 
case, the pairing interaction shifts the chemical potential 
which manifests itself as a time-dependent phase factor 
e -iwt mu itiplying A(i). As discussed in Ref. [5], this can 
be taken into account by a gauge transformation which 
shifts single particle energies by lu/2. In the transformed 
problem, only the modulus |A| varies with time, with its 
functional form still given by the elliptic function (31). 
The parameters A± and ui satisfy algebraic self consis- 
tency equations [5] with tanh ^/3e p replaced by 1 — 2n p . 

The variation in the period T can be linked to the 
fluctuations in the initial state projection on the unsta- 
ble mode (9). Denoting this projection 77, we can write 
for it a distribution of Porter-Thomas form [30] , 



P(\rj\ = x) = 2ux exp(— ux 2 ), 



(44) 



The latter describes fluctuations of individual compo- 
nents of a random complex vector in a high dimensional 
space, with the parameter u being a function of the vec- 
tor norm statistics. In our case, the effective dimension- 
ality can be estimated as a ratio of temperature to the 
level spacing, d ~ l/05e = N/(J3W). To relate rj to the 
period T, we write the time-dependent A at the times 
described by linear instability as A(t) cx rye 7 *. The corre- 
sponding time range can be estimated from the condition 
A(t) < A+ ~ 7, giving t = 7 _1 1117/77. The time t is close 
to the phase lag r which, for the reasons discussed above, 
is approximately equal to ~T. 

Porter-Thomas distribution predicts order of magni- 
tude fluctuations with typical r\ ~ u^ 1 / 2 . This translates 
into fluctuations of T = 2j~ 1 In 7/77 about its mean value 
with the dispersion independent of u. Indeed, Figs. 5, 6 
indicate that a ten-fold increase in temperature, while 
reducing the period, has little effect on its fluctuations. 

To see whether the randomness in the initial state, and 
specifically in the occupancy n p , can explain the fluctu- 
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ations in A + recorded in Figs. 5, 6, we consider instability 
of a particular initial unpaired state (15). Eq.(18) gives 
an equation for the instability exponent: 



9 P 



9 P = 



= |„,(0)|2 _ L,(0)|2 



(45) 



(0) „(o) 



taking values 
= (1-«p) 2 , 



with microscopic non-averaged u p 
zero or one with the probabilities p| Mp |2. 

P\v p \'=i = n\ [see Eq.(15)]. 

The fluctuation Sg p causes deviation <5£ from the av- 
erage value C = *7- Linearization of (45) gives 



Sp 



(2e p - z 7 )2 



E 



2c r 



(46) 



with g p = tanh \f3e P - We see that both the real and the 
imaginary part of <5£ = 5C,' + i8Q" are nonzero, due to the 
parts of the distribution fluctuation Sg p even and odd 
relative to Ep. The imaginary part 8C," gives fluctuation 
in the instability growth exponent 7. The real part 5Q 1 
can be associated with a shift of the chemical potential 
due to particle- hole imbalance in the pair sector. 

We estimate the magnitude of the fluctuations for low 
temperature T <C Aq ~ 7, when Eq.(46) is reduced to 



8Q = -E( 2e P _ *T)^p- 



(47) 



Separating the real and imaginary part, we obtain 

( S (" 2 )=^ 4 4(S 9 2 P ), (SC")=Se^9 2 p )- (48) 
p p 

Here the second moment (<5<7 p ) is given by 



(top) = (K"\ 2 - KT) 2 - (I"p 0) I 2 - K'\ 2 f 

V) 2 + n l - (1 - 2n vf = 2«p(l - »p) 



,(0) | 2 \2 



= (1 - n p ) 2 + n„ - (1 - 2n p ) 2 = 2n p (l - n p ). (49) 



To obtain an order of magnitude estimate, we note that 
the fluctuations Sg p are of order one at |e p | < T and 
exponentially small at |e p | 3> T. Thus we find 

(i) (SC" 2 )^^Se 7 (ii) r>^ft, (50) 
T 

where <5e = W/iV is the level spacing. 

The T 3 ' 2 temperature dependence of the fluctuation 
in 7, Eq.(50.i), can be compared to the distributions of 
the amplitude A + presented in Figs. 5,6. As we demon- 
strated above, A + is numerically close to 7, becoming 
equal to it in the limit of well-separated solitons (see 
Fig. 4). According to the T 3 / 2 law, a ten- fold increase in 
temperature from Ao/100 to Ao/10 should lead to the 
dispersion in 7 increase by a factor of 30, which is indeed 
close to the increase in A + dispersion seen in Figs. 5,6. 
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FIG. 8. Noise suppression at increasing number of states 
TV. The time dependence A(i) recorded from a simulation 
at T = 0.7T C (fi = 2.5A ) for N = 10 2 , 10 3 , 10 4 , 10 5 states, 
with other parameters the same as above. 
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FIG. 9. Temperature dependence of the soliton train am- 
plitude as recorder from the simulation. To suppress noise, 
the number of levels N was gradually increased from N = 10 
at the lowest temperature to TV = 10 5 at the highest tem- 
perature. Analytic fit displays the single 1 / cosh soliton (38) 
amplitude obtained from Eq.(37) at A- = 0. 

The noise, which quickly grows as a function of temper- 
ature, Eq.(50), can be suppressed by reducing the level 
spacing Se = W/N. The dramatic effect of level spacing 
on noise is illustrated in Fig. 8 which presents the time de- 
pendent A(t) at a relatively high temperature T = 0.7T C 
for several values of the number of levels TV. We observe 
that the noise, at this temperature significant at small TV, 
decreases at large TV, with the time dependence assuming 
the elliptic function form (31). As demonstrated in Fig. 9, 
the soliton train amplitude, recorded at TV large enough 
to minimize noise, follows very closely the analytic tem- 
perature dependence, Eq.(37). The value of TV required 
to reduce noise to the level at which the behavior (37) is 



11 



revealed, grows as a function of temperature. While at 
T = as few levels as N = 10 is quite sufficient, we find 
that N increases rapidly as T c is approached. 

To analyze the noise at T close to T c , we again con- 
sider fluctuations in the instability growth rate, given 
by Eq.(46). At these temperatures, since 7 is linear in 
T c — T near T c , we have 7 <C T. In this case, the sum 
E P 9p( 2e P - h)~ 2 in ( 46 ) equals 



4*7 



e p9 P 



((2e p ) 2 +7 2 ) 2 



= ?0 



After averaging over fluctuations Sg p = g p — g p , we ob- 
tain 



) 2 +7 2 ' 



(51) 



Using the expression (49) for (5g p ) at |e p | «T, we have 

/l^|2\ f 8T \\ V- W 16T2 x ,^ 

(Ki) = (^-j ^E (2ep )2 +7 2 = -^-^ ( 52 ) 



regime condition ta <C r e i requires that T c — T is outside 
the region described by the time-dependent Ginzburg- 
Landau equation [3], which is also quite narrow. From 
this one can conclude that the mean field approach, vali- 
dated by (53), remains accurate in an infinite system, at 
least for spatially uniform solutions. 

It is not unconceivable that the mean field theory, 
known to work well in equilibrium at weak coupling, is 
also valid for the dynamical problem with generic spa- 
tially varying A. However, an understanding of this ques- 
tion can only be achieved after the role of spatial fluctu- 
ations in the BCS instability is clarified [22]. Here most 
interesting are spatial fluctuations of the phase of A, and 
the properties of vortices, which presently are not under- 
stood. If the characteristic length scale of phase fluctu- 
ations is of order or larger than the correlation length £, 
which seems likely to be the case, the dynamics of the 
modulus |A| can be obtained in a local approximation, 
using the results of this work and ignoring spatial depen- 
dence. 



V. SPIN 1/2 REPRESENTATION 



By inspecting the right hand side of Eq.(46) we find that 
the fluctuations in 5Q" ', the instability growth rate, dom- 
inate at T close to T c , while the fluctuations in 5Cf , the 
chemical potential, are smaller by a factor 7/T. Thus 
Eq.(52) gives an estimate for the fluctuations in 7 and, 
by the argument used above, also provides an estimate 
of the noise in the soliton train amplitude A + . 

The condition necessary for the noise (52) to be small, 
\5(\ <C 7, translates into 



16T 2 IT 

7T7 3 



(53) 



We see that the minimal level number required to sup- 
press noise grows as (T c — T)~ 3 near the transition. The 
fast growth is consistent with the results of simulation 
presented in Figs. 8,9. 

The level number N, which so far was taken to be arbi- 
trary, can be related to other parameters as follows. For a 
system of size L smaller than the BCS correlation length 
£ = hvp/A, which corresponds to a tightly trapped cold 
gas, N is of the order of the total particle number. This 
can also be written as a relation of N and particle Fermi 
momentum: N ~ (Lpp/2irh) 3 . 

In an infinite system, or in a system of size larger than 
£, one can define an effective N equal to the number of 
particles in the correlation volume, N ~ (£pf /2tt?i) 3 rs 
(Ep/A) 3 . Comparing this to the inequality (53), with 
the identifications W = Ep, 7 ~ A oc T c — T, we obtain 
a condition Ep 3> T 2 , nonrestrictive in the entire inter- 
val < T < T c . Other requirements for the mean field 
approach are also nonrestrictive: (i) The detuning from 
transition, T c — T, must be outside the fluctuation re- 
gion, very narrow at weak coupling; (ii) The collisionless 



In Sec. II we derived Bloch equations (27) from 
Bogoliubov-deGennes equations for u p and v p , by rewrit- 
ing them in the form of Gorkov equations for g p = 
l u p| 2 ~l w p| 2 > /p = 2w p vp, and then recognizing that these 
quantities form a three-component Bloch vector. To gain 
more insight, here we demonstrate a different approach in 
which spin 1/2 operators and Bloch dynamics appear on 
an earlier stage. Following Anderson [31], we define pseu- 
dospins associated with individual Cooper pair states, by 
assigning 'Pauli spin' operators a p = ^(a p ±i<7V) to each 
pair of fcrmion states with opposite momenta as follows 



<T a -Pl 



a -pi a pT 



(54) 



a-p^atp^- This allows to 



and crp = [cr+,CT p ] = a+ T a p 
represent the BCS problem (3) as an ensemble of inter 
acting spins: 



2A£' 



(55) 



where ^ p means a sum over the pairs of states (p, — p). 
Since all the spins interact with each other equally, the 
mean field theory here is exact, just like for the BCS 
problem. The mean field Hamiltonian for each spin is 



T~t P — bp 



a p , b p = (-A',-A",e p ). 



(56) 



Here the z component of the effective field b p , given by 
the single particle energy, is spin-specific, while the trans- 
verse components, the same for all the spins, satisfy 



A^A' + zA" = A^ (<7+>. 

p 



(57) 
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This dynamical selfconsistency relation for time- 
dependent A and <7p~ is identical to the Gorkov equa- 
tion (23). In the ground state each spin is aligned with 
bp, and the spins form a texture near the Fermi surface 
[31], given by Eq.(28), with spin rotation described by 
the Bogoliubov angle. 

The dynamical problem of interest takes the form of a 
Bloch equation for the spins, 

& P = i[U p , <7 P ] = 2b p x cr p (58) 

with the field b p defined selfconsistently by (56), (57). 
Eq. (58), linearized about the texture state, describes col- 
lective excitations of a superconductor with frequency 2A 
[6,31]. Linearized about the unpaired state, Eq. (58) de- 
scribes the BCS instability (10). 

The Hilbert space for the spin Hamiltonian (55) can 
be constructed in a standard fashion, using the states 

W = I- um ■■■> (59) 

as basis vectors, where cr p =t, I correspond to the fully 
occupied and empty pair states. The pair states having 
fermionic occupancy one are to be excluded as they are 
decoupled from the dynamics (55). The Hilbert space 
spanned by the states (59) provides a full representation 
of the Hamiltonian (55). 

The spin states (59), which are identical to the many- 
body pair states (16), provide the most general descrip- 
tion of the problem. Here one can note, however, that 
the mean field relation (57) eliminates dynamical coher- 
ence of different spins. This allows to simplify the state, 
replacing (59) by a product state 

i^=(g)(: p p ). (60) 

Comparing this to the fermionic product states (5), (15), 
we see that the spinor components u p , v p are identical to 
Bogoliubov amplitudes, since the Bogoliubov-dcGcnnes 
dynamics (7) is equivalent to the Bloch dynamics (58). 

One can reduce the spin 1/2 Bloch equations (58) to 
the Bloch equations (27) for classical vectors r p used 
above as follows. Each pair state (p, — p) participates in 
the product (60) with the probability np + ( 1 — n p ) 2 , and 
is excluded with the probability 2n p (l — n p ). Since the 
Bloch equation (58) is linear in <j p , it takes the same form 
when dp is replaced by its expectation value r p = (er p ). 
One can include the probability of having occupancy or 
2 is the expectation value, which makes the norm of r p 
equal np + (1 — n p ) 2 , in agreement with Eq.(26). Thus 
we sec that the spin formulation is indeed equivalent to 
the fermionic formulation employed above. 

VI. DISCUSSION 

This work demonstrates that the unpaired fermionic 
state, after being suddenly presented with pairing inter- 



action, develops a BCS instability which triggers oscilla- 
tions of the pairing amplitude and other quantities. The 
oscillations are periodic in time and are not damped as 
long as particle collisions do not play a role. The oscilla- 
tory behavior comes quite naturally, given that without 
collisions the system cannot lower its energy to that of 
the BCS ground state. 

What comes as a surprise, however, is that the oscil- 
lations have predictable characteristics despite thermal 
noise in the initial conditions. The time-dependent pair- 
ing amplitude is described by the soliton train solutions 
of Jacobi elliptic function dn form [5] in which only the 
parameters such as the period, amplitude and time lag 
depend on the initial conditions. The explanation for 
such a behavior can be traced to the physics of the BCS 
instability. In the latter, when linearization over the un- 
paired state is analyzed, only one mode exhibits instabil- 
ity, while other modes correspond to the perturbations 
that do not grow. As a result, in the evolution of a generic 
unpaired fermion state only the perturbation along the 
unstable direction is amplified by the instability, selecting 
the special elliptic function as a time dependence. The 
accuracy to which the special solution is selected is con- 
troled by the strength of fluctuations in the initial state, 
due to finite temperature and level spacing. 

The selection phenomenon may appear counterintu- 
itive. Here it is instructive to make comparison to the 
results of Ref. [32] which employs the integrability of the 
BCS problem to study time-dependent solutions. The 
large family of solutions obtained in Ref. [32] could leave 
one under impression that they all are equally relevant for 
the evolution of a generic state, which does not agree with 
the results of our simulation and analytic arguments. In- 
stead, as we have seen above, some solutions are singled 
out by the dynamics, while others are not. This peculiar 
situation illustrates that knowing the general solution of 
a nonlinear problem is not necessarily helpful in identify- 
ing the special solution relevant for the physical system 
in which a selection mechanism is at work. 

Let us also comment on some issues of interest not 
considered in this work. One has to do with energy re- 
laxation, left out of the analysis of collisionless BCS dy- 
namics. The relaxation processes relevant for cold gasses 
are due to elastic collisions involving two particle scat- 
tering. The rate of such processes, estimated above as 
Tgj" , is small compared to the typical frequency of os- 
cillations lu ~ A/7i, allowing for undamped oscillations 
over a relatively long time interval < t < r e i. While 
proper treatment of relaxation can only be obtained with 
the help of quantum kinetic theory, one can account for 
it heuristically [5] by inserting a Landau-Lifshitz term in 
the Bloch equation (27), changing b p to 

b p = b p — !p x r P' 

where l p = b p /|b p |. The resulting evolution exhibits 
damped oscillations converging to the BCS ground state 
asymptotically at large times f > t c i. 
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Another phenomenon of interest is related with spa- 
tial fluctuations. Our discussion focused on the dynam- 
ics in a system of finite size, of order or smaller than the 
correlation length, in which we considered spatially uni- 
form A(t), neglecting the spatial dependence altogether. 
In an infinite system, one expects the emergent pairing 
dynamics to exhibit phase fluctuations and vortices si- 
multaneously with the oscillations of the modulus |A|(t) 
similar to that studied above. If the phase fluctuations 
occur at distances larger than the correlation length and 
the vortices are dilute, one can use the estimates made at 
the end of Sec. IV to argue that the modulus |A| will os- 
cillate in a pretty much the same way as in the spatially 
uniform situation. However, more work will be needed 
to clarify this. 



VII. SUMMARY 

In this paper we studied fermionic pairing in a system 
with time-dependent pairing interaction. We analyze the 
situation when the pairing builds up after the interac- 
tion has been abruptly turned on. Theoretical analy- 
sis, supported by numerical simulations, predicts a stage 
of exponential growth, described by BCS instability of 
the unpaired Fermi gas, followed by periodic oscillations 
described by collisionless nonlinear BCS dynamics. We 
consider spatially uniform situation relevant for systems 
of small size, and find that: 

(i) In the collisionless approximation, at times shorter 
than the energy relaxation time, the oscillations are un- 
damped; 

(ii) The time dependence of the pairing amplitude is 
obtained from an exact solution of the nonlinear BCS 
problem [5] , a periodic soliton train described by the Ja- 
cobi elliptic function dn, with parameters depending on 
the microscopic initial conditions; 

(iii) The robustness of the elliptic function behavior is 
explained by a dynamical selection process, in which the 
BCS instability acts to amplify the initial perturbation 
in a specific unstable mode of the system, generating a 
time dependence with predictable characteristics; 

(iv) The fluctuations of the amplitude, period and time 
lag of the soliton train can be accounted for by occupa- 
tion probability fluctuations in the initial state. 
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APPENDIX A 

Here we present yet another derivation of the Bloch 
equation (27), starting from the evolution of time- 
dependent amplitudes u p (t), v p (t). Let us consider the 



Bogoliubov-deGennes equation (7) with the selfconsis- 
tency condition (6). By introducing a new variable 
w p = Up/vp, the pair of linear differential equations (7) 
is reduced to a single nonlinear equation of Riccati form, 



idtWp = 2e p w p + A(t) - A*(i)^ 



(Al) 



which was analyzed in Ref. [19]. The selfconsistency con- 
dition (6), rewritten in terms of w p (t), becomes 



where 



Q P = \u p \ 2 + \v p \ 2 



(A2) 



(A3) 



is the norm of (w p ,i> p ), conserved by the dynamics (7). 
The initial value w p corresponding to the unpaired 

Fermi gas (17) is w p 0) = e^»(l - n p )/n p = e ^ +i ^ , 
with the phases <p p random and uncorrelated for differ- 
ent p. For (u p ,Vp) of the form (17) we have 

Q p = n 2 p + (1 - n p ) 2 

The factor Q p < 1 describes the effect of Pauli blocking 
which was discussed in Sec. II. 




FIG. 10. Stereographic projection (A4) schematic. 

The next step is to perform a reverse stereographic 
projection of the complex variable w, mapping it onto 
the unit sphere r\+r\+ r\ 

2w 



n + ir 2 



1 as follows: 

' 2 -l 



W 



M 2 + i' r3 \W\ 2 + 1 



(A4) 



(see Fig. 10). The dynamics (Al), written in terms of 
z = n + ir 2 and rs, gives 



dz 



— = -2ie p z + 2iAr 3 , = iA* z - iAz* . (A5) 

After rewriting Eq.(A5) in terms of r p = (n, r 2 , r 3 ) p , we 
again obtain the Bloch equation (27), f p = 2b p x r p , 
with the "magnetic field" b p = — (A', A",e p ). 

The selfconsistency condition (A2) takes the form 



p 



(A6) 



The difference in the form of the selfconsistency relations 
(A6 and (23) is due to the difference in normalization of 
r p here and in Sec. II. Here we have |r p | = 1, while in 
Sec. II we had |r p | = np + (1 — n p ) 2 < 1 which is precisely 
the factor Q p needed to account for the difference in the 
norm. 
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APPENDIX B 

Here we estimate the direct heating of the Fermi sys- 
tem due to the shakeup caused by interaction switching 

Wint(i) = A(i) ^ a £ 4 ,a a £ 3 ,/3 a P2 ,/?%>!, « ■ 

Pl+P2=P3+P4 

The interaction time dependence during switching is 
step-like, varying from to A over a characteristic time 
To, so the Fourier component A w = J e tut \(t)dt bah 
Aves as ui^ 1 at o>t <C 1. For example, the inter- 
action switching model \(t) = A(l — e~ t / T °)6(t) gives 
A w = i\j - iujTo)). 

The effective temperature T e ff after switching can be 
estimated from the Golden Rule transition rate for parti- 
cle energy excitation matched by the net energy increase 
and fermion specific heat: 

aT e 2 ff = ^ faniM 1 - n 3)(l -n 4 )|A w | 2 x (Bl) 

u, 1...4 

5(huj - e P3 - e P4 + e Pl + e P2 ) 

with Hi = n Pi the occupation numbers of the states P1...4, 
a = T^v. At T = 0, wc obtain aT c 2 ff = V^ >0 TiwiVjAJ 2 
with = ±j/cj 3 . 

The integral over 10 gives T 2 ff ~ A 2 ^ 3 t ~ 3 . Comparing 
T ff to Ao we find that the system is not overheated, 
T cS < A , when E f t > (An/A ) 2/3 . This condi- 
tion is compatible with the nonadiabaticity requirement 
To <C ta, allowing for a range of possible values of the 
switching times tq for which heating of the Fermi gas is 
negligible. 
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